Relative fitness in a field setting - Figure 4 (Evolution MS)

#Background At Blue Oak Ranch Reserve, we established a 10m x 26m fenced field site to test whether evolution over the course of invasion away from roads has resulted in enhanced performance in undisturbed vegetation relative to roadside populations. The experiment was replicated in a randomized block design (20 plots in total). Each plot was 1.5m2 with 16 D. graveolens growing in a 4 x 4 grid centered on the plot. There was 33cm between each plant and 25cm between the edge plants and the border of the plot.

The experiment included multiple treatments; however, only the two most relevant to the focus of this paper are included here. We tested whether plant genotypes collected from the two habitats (roadside and vegetated) responded differently to the disturbance of biomass removal. We tilled in December 2020 to completely remove below and aboveground biomass, and then weeded to remove aboveground biomass throughout the growing season. In contrast, we left the control plots untouched, allowing the previous year’s thatch to persist and background vegetation to grow throughout the experiment.

In January 2021, we germinated seeds in Petri dishes at the UCSC Coastal Science Campus greenhouse incubation chambers before transplanting them into field-collected soil (collected in late December 2020 from Blue Oak Ranch Reserve). Seedlings grew in the greenhouse for about eight weeks until all plants had their first two true leaves emerge and lengthen. Ideally, we would have placed seeds directly into the field, but to maximize biosafety, we used seedling transplants that could be tracked with 100% certainty.

We measured the longest leaf for each plant and then transplanted them into the ground in late February 2021 at Blue Oak Ranch Reserve. During the first month of growth, we replaced any D. graveolens that died. We conducted weekly phenology surveys to assess D. graveolens plant health, and at the first sign of buds, we measured plant height and harvested the aboveground biomass by cutting at the root crown and drying in a 60ºC oven for 3 days before weighing.

#Data Analysis Statistical analyses were performed in R version 4.2.1 (R Core Team 2022) using linear mixed-effects models with the lme4 (Bates et al. 2015), lmerTest (Kuznetsova et al. 2017), and DHARMa packages (Hartig 2022), generalized linear mixed models with the glmmTMB package (Brooks et al. 2017), and mixed effects cox models with the coxme (Therneau 2022a) and survival (Therneau 2022b) packages.

#Libraries

#install.packages("coxme")
#install.packages("survival")
#install.packages("ggplot2")
#install.packages("ggfortify")
#install.packages("car")
#install.packages("multcomp")
#install.packages("lme4")
#install.packages("lmerTest")
#install.packages("DHARMa")
#install.packages("dplyr")
#install.packages("emmeans")
#install.packages('TMB', type = 'source')
#install.packages("glmmTMB")
#install.packages("MASS")
#install.packages("emmeans")
#install.packages("AICcmodavg")
library(coxme)
library(survival)
library(ggplot2)
library(ggfortify)
library(car)
library(multcomp)
library(lme4)
library(lmerTest)
library(DHARMa)
library(dplyr)
library(emmeans)
library(TMB)
library(glmmTMB)
library(MASS)
library(emmeans)
library(AICcmodavg)
library(tidyverse)

#Load Data This dataframe has one row per plant (320 observations). Data are for survivorship curves (3 censor options), the number of days the plant stayed alive (NumDaysAlive) and aboveground biomass. Censors with a 1 denote reaching the event (CensorAll = died, CensorBiomass = survived to collect biomass, CensorReproduction = survived to reproduce) and a 0 denoting when a seed didn’t germinate by the last census date (Census = 11/15/21). CensorReproduction will be most useful in understanding the amount of biomass produced by an individual when buds appear.

mydata<-read.csv("/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/data/field_relative_fitness/MS_Evolution_Fig4_2021.csv",stringsAsFactors=T)
str(mydata) #Check that each column has the right class (factor, integer, numeric, etc.)

'data.frame':   320 obs. of  35 variables:
 $ Block             : Factor w/ 10 levels "A","B","C","D",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ Plot              : int  5 2 2 1 5 1 1 4 4 2 ...
 $ Flag              : Factor w/ 20 levels "A2","A5","B1",..: 2 4 5 7 10 11 13 16 17 19 ...
 $ Pos               : int  3 15 13 13 6 10 8 4 12 3 ...
 $ Flag_Pos          : Factor w/ 320 levels "A2_01","A2_02",..: 19 63 77 109 150 170 200 244 268 291 ...
 $ Population        : Factor w/ 16 levels "BAY-A","BAY-O",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Site              : Factor w/ 8 levels "BAY","CHE","GUA",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Habitat           : Factor w/ 2 levels "Roadside","Vegetated": 1 1 1 1 1 1 1 1 1 1 ...
 $ Treatment         : Factor w/ 2 levels "Biomass Removal",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ HabitatTreatment  : Factor w/ 4 levels "Roadside/Biomass Removal",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ PlantDate         : Factor w/ 7 levels "2/27/21","3/1/21",..: 5 1 1 1 1 1 7 1 1 1 ...
 $ LeafMeas1         : num  30 20.3 27 14.9 16.7 ...
 $ LeafMeas2         : int  46 51 NA 56 42 45 36 58 52 45 ...
 $ Growth            : num  16 30.7 -27 41 25.3 ...
 $ MortDate          : Factor w/ 27 levels "","10/1/21","11/14/21",..: 1 24 7 1 1 1 1 1 1 1 ...
 $ MortHarvDate      : Factor w/ 32 levels "10/1/21","10/8/21",..: 20 26 8 28 19 28 28 28 29 29 ...
 $ CensorAll         : int  1 1 1 1 1 1 1 1 1 1 ...
 $ DaysMort          : int  105 181 62 195 139 195 185 195 198 198 ...
 $ Census            : Factor w/ 1 level "11/15/21": 1 1 1 1 1 1 1 1 1 1 ...
 $ DaysCensus        : int  241 261 261 261 261 261 251 261 261 261 ...
 $ NumDaysAlive      : int  105 181 62 195 139 195 185 195 198 198 ...
 $ HarvDate          : Factor w/ 19 levels "","10/1/21","10/8/21",..: 7 13 1 15 6 15 15 15 16 16 ...
 $ CensorBiomass     : int  1 1 0 1 1 1 1 1 1 1 ...
 $ BudDate           : Factor w/ 18 levels "","10/1/21","10/8/21",..: 7 1 1 15 9 15 15 15 15 15 ...
 $ CensorReproduction: int  1 0 0 1 1 1 1 1 1 1 ...
 $ SurvToRepro       : int  1 0 0 1 1 1 1 1 1 1 ...
 $ CensorTest        : int  0 1 1 0 0 0 0 0 0 0 ...
 $ PropBud           : num  0.763 0.763 0.763 0.763 0.763 0.763 0.763 0.763 0.763 0.763 ...
 $ PropBudSite       : num  0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 ...
 $ Phenology         : Factor w/ 31 levels "10/1/21","10/8/21",..: 17 26 8 28 19 28 28 28 28 28 ...
 $ DaysToPheno       : int  98 181 62 195 139 195 185 195 195 195 ...
 $ CensorPheno       : int  0 1 1 0 0 0 0 0 0 0 ...
 $ Height            : num  12.4 33.8 NA 56.6 21.5 55.2 45.2 37.5 50.8 48.5 ...
 $ Biomass           : num  0.762 7.048 NA 7.803 2.572 ...
 $ Biomass.date      : Factor w/ 16 levels "","10/13/21",..: 15 14 1 6 13 9 9 9 8 8 ...

mydata$Site<-as.character(mydata$Site)
mydata$Treatment<-factor(mydata$Treatment,
                         levels=c("Grassland","Biomass Removal")) #Changing the contrast order so that everything is compared to Grassland (control)

#Early Growth This code uses Growth data with Habitat (roadside and vegetated) and Treatment in a glmm model. Anova and Tukey tests are used on the successful Model4 with the creation of a box plot as a finished product.

Note: 10 blocks, 2 treatments, 16 populations (CHE-O), 8 sites (random: population pairs; CHE), 2 habitat (fixed effect: roadside and vegetated; R or O)

Also Note: Growth data is a measure of growth of longest leaf. Each plant was measured upon transplanting into the field in March and then again in May 2021. This plant has a juvenile stage of a basal rosette, and then it bolts and produced smaller cauline leaves. The goal was to capture early growth data for this plant before bolting occurs so that the measurements capture the basal rosette stage, however in some cases the plants bolted earlier than expected and the resulting measurement was smaller than previous. In these cases we determined that the data should be removed from the dataset as the negative number (or changing it to a zero) does not reflect the biological importance of the measurement.

Here we need to filter the data to remove Growth>0 and to convert plant measurement dates to date format

growth_mydata<-mydata%>%filter(Growth>0) #Here I am only looking at the Growth data that is greater than 0 (see Also Note above)
growth_mydata$PlantDate<-as.Date(growth_mydata$PlantDate,"%m/%d/%y") 

growth_mydata$Num.Days.Growth<-as.Date("2021-05-22")-growth_mydata$PlantDate #number days
growth_mydata$Num.Days.Growth<-as.numeric(growth_mydata$Num.Days.Growth)
growth_mydata$Growth.Rate<-growth_mydata$Growth/growth_mydata$Num.Days.Growth #First calculate number of days of growth to get the Rate 

##Histograms ###Original data When we plot the original data as a histogram, we find that it is skewed. We also use the Shapiro-Wilk normality test, but get a p-value = 1.402e-05.

In statistical hypothesis testing, a common significance level (alpha) is set at 0.05. If the p-value is less than alpha (p-value < alpha), it suggests strong evidence against the null hypothesis. In this case, the p-value is extremely small (less than 1.402e-05), indicating very strong evidence against the null hypothesis.

Therefore, based on the given output, it can be concluded that the dataset did not pass the Shapiro-Wilk normality test. The data is unlikely to follow a normal distribution. Let’s consider log transforming the data.

hist(growth_mydata$Growth.Rate,
     col='steelblue',
     main='Original') 

shapiro.test(growth_mydata$Growth.Rate)

    Shapiro-Wilk normality test

data:  growth_mydata$Growth.Rate
W = 0.96321, p-value = 1.402e-05

###Log transform data (https://www.statology.org/transform-data-in-r/) When we log transform the data and plot using a histogram, the data does not look, but we need to still test for normalicy with a Shapiro-Wilk normality test. Here we find that the data are not normally distributed because the p-value = 7.398e-14.

log_growth_mydata<-log10(growth_mydata$Growth.Rate)
hist(log_growth_mydata,
     col='steelblue',main='Log Transformed') #Log transformed data does not look good at all!

shapiro.test(log_growth_mydata) #Data does not improve with log transformation

    Shapiro-Wilk normality test

data:  log_growth_mydata
W = 0.85341, p-value = 7.398e-14

qqnorm(log_growth_mydata) #qqplot
qqline(log_growth_mydata) #add the line 

###Square root transform data This improved the distribution, but it failed the Shapiro-Wilk normality test (p-value = 0.001034). Let’s try poisson.

sqrt_growth_mydata<-sqrt(growth_mydata$Growth.Rate)
hist(sqrt_growth_mydata,
     col='steelblue',main='Square Root Transformed') #Square root transformed data looks better than the original distribution

shapiro.test(sqrt_growth_mydata)

    Shapiro-Wilk normality test

data:  sqrt_growth_mydata
W = 0.97721, p-value = 0.001034

qqnorm(sqrt_growth_mydata) #qqplot
qqline(sqrt_growth_mydata) #add the line 

###Poisson distribution Poisson is actually not a good fit because the growth rate is not an integer (has decimals), which poisson is not equipt to deal with. Maybe I’ll try gamma?

library(MASS)
MASS::fitdistr((growth_mydata$Growth.Rate*1000),
               "Poisson")

Warning: non-integer x = 365.714286Warning: non-integer x = 488.690476Warning: non-integer x = 301.547619Warning: non-integer x = 224.523810Warning: non-integer x = 216.216216Warning: non-integer x = 486.666667Warning: non-integer x = 350.119048Warning: non-integer x = 393.928571Warning: non-integer x = 213.333333Warning: non-integer x = 352.619048Warning: non-integer x = 299.642857Warning: non-integer x = 329.268293Warning: non-integer x = 317.261905Warning: non-integer x = 343.571429Warning: non-integer x = 469.642857Warning: non-integer x = 270.270270Warning: non-integer x = 381.190476Warning: non-integer x = 84.745763Warning: non-integer x = 176.428571Warning: non-integer x = 238.809524Warning: non-integer x = 152.619048Warning: non-integer x = 440.677966Warning: non-integer x = 472.380952Warning: non-integer x = 452.500000Warning: non-integer x = 443.452381Warning: non-integer x = 90.238095Warning: non-integer x = 384.761905Warning: non-integer x = 101.666667Warning: non-integer x = 588.571429Warning: non-integer x = 558.809524Warning: non-integer x = 133.571429Warning: non-integer x = 253.521127Warning: non-integer x = 295.774648Warning: non-integer x = 511.666667Warning: non-integer x = 529.404762Warning: non-integer x = 70.422535Warning: non-integer x = 137.976190Warning: non-integer x = 171.875000Warning: non-integer x = 220.476190Warning: non-integer x = 288.135593Warning: non-integer x = 253.731343Warning: non-integer x = 456.071429Warning: non-integer x = 438.095238Warning: non-integer x = 576.271186Warning: non-integer x = 178.333333Warning: non-integer x = 524.285714Warning: non-integer x = 251.309524Warning: non-integer x = 206.904762Warning: non-integer x = 472.972973Warning: non-integer x = 496.428571Warning: non-integer x = 512.619048Warning: non-integer x = 308.690476Warning: non-integer x = 432.023810Warning: non-integer x = 433.809524Warning: non-integer x = 384.761905Warning: non-integer x = 215.952381Warning: non-integer x = 403.452381Warning: non-integer x = 267.738095Warning: non-integer x = 321.309524Warning: non-integer x = 585.365854Warning: non-integer x = 373.452381Warning: non-integer x = 339.166667Warning: non-integer x = 205.476190Warning: non-integer x = 263.928571Warning: non-integer x = 687.261905Warning: non-integer x = 314.166667Warning: non-integer x = 207.619048Warning: non-integer x = 321.071429Warning: non-integer x = 571.547619Warning: non-integer x = 194.523810Warning: non-integer x = 148.928571Warning: non-integer x = 372.738095Warning: non-integer x = 459.285714Warning: non-integer x = 394.523810Warning: non-integer x = 352.112676Warning: non-integer x = 265.595238Warning: non-integer x = 394.642857Warning: non-integer x = 722.500000Warning: non-integer x = 189.642857Warning: non-integer x = 327.261905Warning: non-integer x = 546.875000Warning: non-integer x = 498.452381Warning: non-integer x = 448.690476Warning: non-integer x = 337.837838Warning: non-integer x = 243.243243Warning: non-integer x = 157.500000Warning: non-integer x = 539.047619Warning: non-integer x = 129.880952Warning: non-integer x = 427.619048Warning: non-integer x = 546.875000Warning: non-integer x = 448.928571Warning: non-integer x = 313.690476Warning: non-integer x = 328.358209Warning: non-integer x = 119.047619Warning: non-integer x = 530.357143Warning: non-integer x = 491.666667Warning: non-integer x = 432.835821Warning: non-integer x = 474.576271Warning: non-integer x = 578.125000Warning: non-integer x = 247.023810Warning: non-integer x = 221.666667Warning: non-integer x = 291.428571Warning: non-integer x = 310.810811Warning: non-integer x = 131.071429Warning: non-integer x = 435.833333Warning: non-integer x = 498.809524Warning: non-integer x = 278.095238Warning: non-integer x = 250.357143Warning: non-integer x = 499.642857Warning: non-integer x = 204.047619Warning: non-integer x = 350.595238Warning: non-integer x = 329.285714Warning: non-integer x = 320.357143Warning: non-integer x = 390.357143Warning: non-integer x = 379.404762Warning: non-integer x = 396.071429Warning: non-integer x = 237.261905Warning: non-integer x = 582.738095Warning: non-integer x = 132.976190Warning: non-integer x = 420.119048Warning: non-integer x = 432.835821Warning: non-integer x = 475.609756Warning: non-integer x = 268.928571Warning: non-integer x = 231.707317Warning: non-integer x = 81.547619Warning: non-integer x = 546.875000Warning: non-integer x = 253.809524Warning: non-integer x = 357.619048Warning: non-integer x = 362.857143Warning: non-integer x = 141.785714Warning: non-integer x = 64.880952Warning: non-integer x = 610.169492Warning: non-integer x = 322.380952Warning: non-integer x = 432.619048Warning: non-integer x = 343.809524Warning: non-integer x = 324.324324Warning: non-integer x = 310.810811Warning: non-integer x = 457.142857Warning: non-integer x = 538.452381Warning: non-integer x = 407.261905Warning: non-integer x = 397.976190Warning: non-integer x = 121.621622Warning: non-integer x = 81.081081Warning: non-integer x = 59.701493Warning: non-integer x = 140.357143Warning: non-integer x = 95.238095Warning: non-integer x = 21.785714Warning: non-integer x = 17.738095Warning: non-integer x = 189.642857Warning: non-integer x = 104.477612Warning: non-integer x = 4.880952Warning: non-integer x = 173.333333Warning: non-integer x = 44.642857Warning: non-integer x = 1.785714Warning: non-integer x = 168.214286Warning: non-integer x = 100.714286Warning: non-integer x = 16.190476Warning: non-integer x = 67.567568Warning: non-integer x = 74.626866Warning: non-integer x = 57.500000Warning: non-integer x = 243.333333Warning: non-integer x = 282.023810Warning: non-integer x = 193.809524Warning: non-integer x = 31.190476Warning: non-integer x = 140.357143Warning: non-integer x = 8.214286Warning: non-integer x = 214.404762Warning: non-integer x = 91.904762Warning: non-integer x = 4.285714Warning: non-integer x = 23.214286Warning: non-integer x = 68.809524Warning: non-integer x = 60.476190Warning: non-integer x = 147.023810Warning: non-integer x = 135.357143Warning: non-integer x = 52.142857Warning: non-integer x = 143.214286Warning: non-integer x = 29.850746Warning: non-integer x = 195.476190Warning: non-integer x = 2.380952Warning: non-integer x = 89.552239Warning: non-integer x = 169.047619Warning: non-integer x = 60.119048Warning: non-integer x = 194.642857Warning: non-integer x = 57.023810Warning: non-integer x = 3.809524Warning: non-integer x = 48.452381Warning: non-integer x = 26.309524Warning: non-integer x = 189.189189Warning: non-integer x = 93.571429Warning: non-integer x = 136.428571Warning: non-integer x = 21.904762Warning: non-integer x = 33.928571Warning: non-integer x = 168.333333Warning: non-integer x = 152.542373Warning: non-integer x = 240.238095Warning: non-integer x = 134.047619Warning: non-integer x = 175.675676Warning: non-integer x = 176.309524Warning: non-integer x = 247.142857Warning: non-integer x = 66.666667Warning: non-integer x = 68.333333Warning: non-integer x = 243.243243Warning: non-integer x = 109.756098Warning: non-integer x = 72.976190Warning: non-integer x = 43.928571Warning: non-integer x = 16.904762Warning: non-integer x = 304.523810Warning: non-integer x = 87.142857Warning: non-integer x = 64.880952Warning: non-integer x = 76.071429Warning: non-integer x = 189.189189Warning: non-integer x = 54.054054Warning: non-integer x = 44.776119Warning: non-integer x = 98.214286Warning: non-integer x = 219.512195Warning: non-integer x = 87.142857Warning: non-integer x = 9.642857Warning: non-integer x = 117.976190Warning: non-integer x = 41.071429Warning: non-integer x = 214.166667

     lambda  
  264.890459 
 (  1.082627)

qqPlot(growth_mydata$Growth.Rate,
       distribution="pois",
       lambda=1)

[1] 40 81

###Gamma distribution Nope, not gamma

gamma<-fitdistr(growth_mydata$Growth.Rate,
                "gamma")

Warning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs produced

qqp(growth_mydata$Growth.Rate,
    "gamma",
    shape=gamma$estimate[[1]],
    rate=gamma$estimate[[2]])

[1] 40 81

##Models ###Full Model 1 - glmm

growth_fullmodel1<-glmmTMB(Growth.Rate~
                             Habitat*Treatment+
                             (1|Site)+
                             (1|Block),
                           family=Gamma(link="log"),
                           data=growth_mydata) 
summary(growth_fullmodel1)

 Family: Gamma  ( log )
Formula:          Growth.Rate ~ Habitat * Treatment + (1 | Site) + (1 | Block)
Data: growth_mydata

     AIC      BIC   logLik deviance df.resid 
  -285.7   -261.8    149.9   -299.7      219 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev.
 Site   (Intercept) 0.0008935 0.02989 
 Block  (Intercept) 0.0116777 0.10806 
Number of obs: 226, groups:  Site, 8; Block, 10

Dispersion estimate for Gamma family (sigma^2): 0.392 

Conditional model:
                                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)                               -2.38870    0.11451 -20.860   <2e-16 ***
HabitatVegetated                           0.14981    0.14068   1.065    0.287    
TreatmentBiomass Removal                   1.32254    0.13255   9.977   <2e-16 ***
HabitatVegetated:TreatmentBiomass Removal -0.09731    0.17488  -0.556    0.578    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

qqnorm(resid(growth_fullmodel1)) #qqplot
qqline(resid(growth_fullmodel1)) #add the line

testDispersion(growth_fullmodel1) #red line should be in the middle of the distribution

    DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated

data:  simulationOutput
dispersion = 0.434, p-value < 2.2e-16
alternative hypothesis: two.sided

myDHARMagraph_grow1<-simulateResiduals(growth_fullmodel1) #making a graph using DHARMa package, also testing for heteroscedasticity, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html#heteroscedasticity
plot(myDHARMagraph_grow1) #plotting graph. Looks good, but check the outliers to make sure they are real.

###Full Model 2 - lmer

growth_fullmodel2<-lmer(sqrt(Growth.Rate)~
                          Habitat*Treatment+
                          (1|Site)+
                          (1|Block),
                        data=growth_mydata)
isSingular(growth_fullmodel2,
           tol=1e-4) #=False

[1] FALSE

summary(growth_fullmodel2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: sqrt(Growth.Rate) ~ Habitat * Treatment + (1 | Site) + (1 | Block)
   Data: growth_mydata

REML criterion at convergence: -271.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.75170 -0.70742  0.04675  0.72993  2.37882 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Block    (Intercept) 0.0009601 0.03099 
 Site     (Intercept) 0.0001233 0.01110 
 Residual             0.0154163 0.12416 
Number of obs: 226, groups:  Block, 10; Site, 8

Fixed effects:
                                            Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)                                 0.274343   0.022973  60.272607   11.94   <2e-16 ***
HabitatVegetated                            0.026520   0.027903 210.393040    0.95    0.343    
TreatmentBiomass Removal                    0.297160   0.025038 214.982016   11.87   <2e-16 ***
HabitatVegetated:TreatmentBiomass Removal  -0.008642   0.034637 209.501173   -0.25    0.803    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) HbttVg TrtmBR
HabittVgttd -0.634              
TrtmntBmssR -0.723  0.581       
HbttVgt:TBR  0.509 -0.804 -0.711

anova(growth_fullmodel2)

Type III Analysis of Variance Table with Satterthwaite's method
                  Sum Sq Mean Sq NumDF  DenDF  F value Pr(>F)    
Habitat           0.0252  0.0252     1 210.09   1.6364 0.2022    
Treatment         4.2607  4.2607     1 220.63 276.3728 <2e-16 ***
Habitat:Treatment 0.0010  0.0010     1 209.50   0.0623 0.8032    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

qqnorm(resid(growth_fullmodel2)) #qqplot
qqline(resid(growth_fullmodel2)) #add the line

testDispersion(growth_fullmodel2) #red line should be in the middle of the distribution

    DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated

data:  simulationOutput
dispersion = 0.98259, p-value = 0.92
alternative hypothesis: two.sided

myDHARMagraph_grow2<-simulateResiduals(growth_fullmodel2) #making a graph using DHARMa package, also testing for heteroscedasticity, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html#heteroscedasticity
plot(myDHARMagraph_grow2) #plotting graph. At this point, you don't want any text or lines to be red. The QQ plots are not looking good for this model. Let's try fitting to a negative binomial distribution.

###Other Model Attempts ####lmer with log data Modeling Habitat and Treatment with Site and Block as random effects. For this model, we use the log(Growth.Rate) data. The qqplot does not look good.

####glmer.nb: Negative Binomial Here I tried fitting my data to a negative binomial distribution, but my data are non-integer, so this model won’t work. Let’s try building a model and fitting it to a Beta distribution which can deal with non-integers (between 0 and 1), first without a link function. Check it with DHARMa, and if it doesn’t look good, then fit it to a Beta distribution with a logit link function.

fullmodel<-glmer.nb(Growth.Rate)~ #Response variable: growth rate HabitatTreatment+ #Fixed effects and their interactions() (1|Site)+(1|Block), #Random effect with random intercept only data=mydata) #Dataframe

Generalized linear mixed-effects model (GLMM) with a negative binomial distribution is part of the lme4 package, which is commonly used for fitting various types of mixed-effects models.

The negative binomial distribution is often used to model count data that exhibit overdispersion, meaning the variance is greater than the mean. This distribution is useful when the assumption of a Poisson distribution (equal mean and variance) is violated.

The glmer.nb() function specifically fits a GLMM with a negative binomial distribution by taking into account both fixed effects (predictors) and random effects (grouping factors). It allows for the modeling of correlated data where observations within the same group may be more similar to each other than to observations from other groups.

####glmm: Beta Distribution Turns out this was also not helpful because beta distributions are for proportions, which this is not.But I’ll keep this here for future reference.

fullmodel<-glmmTMB(Growth.Rate~ #Response variable: growth rate HabitatTreatment+ #Fixed effects and their interactions () (1|Site)+(1|Block), #Random effect with random intercept only family=beta_family(), #Specify the family as beta_family (for beta regression) data=growth_mydata) #Dataframe

Generalized linear mixed-effects models (GLMMs) using a template model builder (TMB) framework. GLMMs are a type of statistical model that extends generalized linear models (GLMs) to account for correlated or clustered data, hierarchical structures, and random effects. They are suitable for analyzing data with non-normal response variables or data that violate the assumptions of traditional linear models.

The glmmTMB() function allows you to specify a GLMM using a formula syntax similar to GLMs. It supports various families of distributions (e.g., Gaussian, binomial, Poisson) and link functions, allowing for the analysis of different types of response variables. Additionally, it allows for the inclusion of both fixed effects (predictors) and random effects (grouping factors) to capture the variability within and between groups.

The TMB framework implemented in glmmTMB utilizes efficient algorithms for model fitting and parameter estimation. It offers computational advantages, making it particularly useful for large datasets or models with complex structures. The package also provides functionality for handling zero-inflation, overdispersion, and handling non-Gaussian response distributions.

###Best Model

growth_fullmodel2<-lmer(sqrt(Growth.Rate)~
                          Habitat*Treatment+
                          (1|Site)+
                          (1|Block),
                        data=growth_mydata)
summary(growth_fullmodel2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: sqrt(Growth.Rate) ~ Habitat * Treatment + (1 | Site) + (1 | Block)
   Data: growth_mydata

REML criterion at convergence: -271.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.75170 -0.70742  0.04675  0.72993  2.37882 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Block    (Intercept) 0.0009601 0.03099 
 Site     (Intercept) 0.0001233 0.01110 
 Residual             0.0154163 0.12416 
Number of obs: 226, groups:  Block, 10; Site, 8

Fixed effects:
                                            Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)                                 0.274343   0.022973  60.272607   11.94   <2e-16 ***
HabitatVegetated                            0.026520   0.027903 210.393040    0.95    0.343    
TreatmentBiomass Removal                    0.297160   0.025038 214.982016   11.87   <2e-16 ***
HabitatVegetated:TreatmentBiomass Removal  -0.008642   0.034637 209.501173   -0.25    0.803    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) HbttVg TrtmBR
HabittVgttd -0.634              
TrtmntBmssR -0.723  0.581       
HbttVgt:TBR  0.509 -0.804 -0.711

#Survival Analysis This code uses NumDaysAlive data with Habitat (roadside and vegetated) and Treatment in a Cox proportional hazards model to assess survival. Information here: https://cran.r-project.org/web/packages/coxme/vignettes/coxme.pdf. “Surv” creates a survival object to combine the days column (NumDaysAlive) and the reproductive censor column (ReproductionCensor) to be used as a response variable in a model formula. The model follows the same syntax as linear models (lm, lmer, etc). Fixed effects: “var1 * var2” will give you the interaction term and the individual variables: so “var1 + var2 + var1:var2”. To add random effects, type “+ (1|random effect)”.

##Histograms

#Number of Days Alive - All data
hist(mydata$NumDaysAlive,
     col='steelblue',
     main='Original') 

#Number of Days Alive - By Habitat
ggplot(mydata,
       aes(x=NumDaysAlive))+
  geom_histogram()+
  facet_wrap(vars(Habitat)) #Here we see that both Habitats have the same bi-modal distribution.

#Number of Days Alive - By Treatment
ggplot(mydata,
       aes(x=NumDaysAlive))+
  geom_histogram()+
  facet_wrap(vars(Treatment))

##Models ###Full Model 1 - glm - generalized linear model

###Cox Models Cox proportional hazard models Here we will use ‘coxme’ which allows you to conduct mixed effects Cox proportional hazards models. Information here: https://cran.r-project.org/web/packages/coxme/vignettes/coxme.pdf. We conducted a germination experiment using Dittrichia graveolens seeds on filter paper. “Surv” creates a survival object to combine the days column (NumDaysAlive) and the censor column (CensorReproduction) to be used as a response variable in a model formula. The model follows the same syntax as linear models (lm, lmer, etc). Fixed effects: “var1 * var2” will give you the interaction term and the individual variables: so “var1 + var2 + var1:var2”. To add random effects, type “+ (1|random effect)”.

Assumptions for cox models: https://www.theanalysisfactor.com/assumptions-cox-regression/

####Simple Model Start by making a simple model with no random effects. This will be compared to the full model with random effects.

cox_simplemodel1<-coxph(Surv(NumDaysAlive,
                             CensorReproduction)~ #Only those that survived to bud are included
                          Habitat*Treatment,
                        data=mydata)
summary(cox_simplemodel1)

Call:
coxph(formula = Surv(NumDaysAlive, CensorReproduction) ~ Habitat * 
    Treatment, data = mydata)

  n= 320, number of events= 160 

                                             coef exp(coef) se(coef)      z Pr(>|z|)    
HabitatVegetated                          -0.5771    0.5615   0.3317 -1.740   0.0819 .  
TreatmentBiomass Removal                   1.6992    5.4694   0.3472  4.894 9.89e-07 ***
HabitatVegetated:TreatmentBiomass Removal  0.5756    1.7782   0.3781  1.523   0.1279    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                                          exp(coef) exp(-coef) lower .95 upper .95
HabitatVegetated                             0.5615     1.7809    0.2931     1.076
TreatmentBiomass Removal                     5.4694     0.1828    2.7694    10.802
HabitatVegetated:TreatmentBiomass Removal    1.7782     0.5624    0.8476     3.731

Concordance= 0.586  (se = 0.035 )
Likelihood ratio test= 67.19  on 3 df,   p=2e-14
Wald test            = 43.68  on 3 df,   p=2e-09
Score (logrank) test = 53.87  on 3 df,   p=1e-11

print(cox_simplemodel1)

Call:
coxph(formula = Surv(NumDaysAlive, CensorReproduction) ~ Habitat * 
    Treatment, data = mydata)

                                             coef exp(coef) se(coef)      z        p
HabitatVegetated                          -0.5771    0.5615   0.3317 -1.740   0.0819
TreatmentBiomass Removal                   1.6992    5.4694   0.3472  4.894 9.89e-07
HabitatVegetated:TreatmentBiomass Removal  0.5756    1.7782   0.3781  1.523   0.1279

Likelihood ratio test=67.19  on 3 df, p=1.707e-14
n= 320, number of events= 160 

predict(cox_simplemodel1)

  [1]  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676
  [9]  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676
 [17]  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676
 [25]  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676
 [33]  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676
 [41]  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676
 [49]  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676
 [57]  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676
 [65]  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676
 [73]  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676  1.6991676
 [81]  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849
 [89]  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849
 [97]  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849
[105]  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849
[113]  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849
[121]  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849
[129]  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849
[137]  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849
[145]  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849
[153]  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849  1.6976849
[161]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
[169]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
[177]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
[185]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
[193]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
[201]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
[209]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
[217]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
[225]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
[233]  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000
[241] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[249] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[257] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[265] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[273] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[281] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[289] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[297] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[305] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932
[313] -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932 -0.5770932

hist(predict(cox_simplemodel1))

####Full Model 1 - coxme Now make a full model using random effects

cox_fullmodel1<-coxme(Surv(NumDaysAlive,CensorReproduction)~
                        Habitat*Treatment+
                        (1|Site)+
                        (1|Block),
                      data=mydata)
summary(cox_fullmodel1)

Cox mixed-effects model fit by maximum likelihood
  Data: mydata
  events, n = 160, 320
  Iterations= 7 47 
                    NULL Integrated    Fitted
Log-likelihood -676.4353  -623.8723 -606.6463

                   Chisq   df p    AIC   BIC
Integrated loglik 105.13  5.0 0  95.13 79.75
 Penalized loglik 139.58 13.3 0 112.98 72.07

Model:  Surv(NumDaysAlive, CensorReproduction) ~ Habitat * Treatment +      (1 | Site) + (1 | Block) 
Fixed coefficients
                                                coef exp(coef)  se(coef)     z       p
HabitatVegetated                          -0.4439054 0.6415261 0.3492953 -1.27 2.0e-01
TreatmentBiomass Removal                   1.9078245 6.7384132 0.3776068  5.05 4.4e-07
HabitatVegetated:TreatmentBiomass Removal  0.5337000 1.7052300 0.3927154  1.36 1.7e-01

Random effects
 Group Variable  Std Dev    Variance  
 Site  Intercept 0.15813763 0.02500751
 Block Intercept 0.81215553 0.65959660

print(cox_fullmodel1)

Cox mixed-effects model fit by maximum likelihood
  Data: mydata
  events, n = 160, 320
  Iterations= 7 47 
                    NULL Integrated    Fitted
Log-likelihood -676.4353  -623.8723 -606.6463

                   Chisq   df p    AIC   BIC
Integrated loglik 105.13  5.0 0  95.13 79.75
 Penalized loglik 139.58 13.3 0 112.98 72.07

Model:  Surv(NumDaysAlive, CensorReproduction) ~ Habitat * Treatment +      (1 | Site) + (1 | Block) 
Fixed coefficients
                                                coef exp(coef)  se(coef)     z       p
HabitatVegetated                          -0.4439054 0.6415261 0.3492953 -1.27 2.0e-01
TreatmentBiomass Removal                   1.9078245 6.7384132 0.3776068  5.05 4.4e-07
HabitatVegetated:TreatmentBiomass Removal  0.5337000 1.7052300 0.3927154  1.36 1.7e-01

Random effects
 Group Variable  Std Dev    Variance  
 Site  Intercept 0.15813763 0.02500751
 Block Intercept 0.81215553 0.65959660

predict(cox_fullmodel1)

  [1]  3.570525077  1.861119361  1.651958824  1.760682078  3.024951374  2.034175698  2.020368065
  [8]  1.892909583  1.145348506  1.050083552  3.525864067  1.816458351  1.607297813  1.716021068
 [15]  2.980290364  1.989514688  1.975707055  1.848248573  1.100687496  1.005422542  3.405466903
 [22]  1.696061187  1.486900650  1.595623904  2.859893200  1.869117524  1.855309892  1.727851409
 [29]  0.980290332  0.885025378  3.506920692  1.797514976  1.588354439  1.697077693  2.961346989
 [36]  1.970571313  1.956763680  1.829305198  1.081744121  0.986479167  3.424679661  1.715273945
 [43]  1.506113408  1.614836662  2.879105958  1.888330282  1.874522649  1.747064167  0.999503090
 [50]  0.904238136  3.610092277  1.900686561  1.691526024  1.800249278  3.064518574  2.073742898
 [57]  2.059935265  1.932476783  1.184915706  1.089650752  3.363249770  1.653844054  1.444683517
 [64]  1.553406771  2.817676067  1.826900391  1.813092758  1.685634276  0.938073199  0.842808245
 [71]  3.410300187  1.700894471  1.491733934  1.600457188  2.864726484  1.873950808  1.860143175
 [78]  1.732684693  0.985123616  0.889858662  3.660319680  1.950913964  1.741753427  1.850476681
 [85]  3.114745977  2.123970301  2.110162668  1.982704186  1.235143109  1.139878155  3.615658670
 [92]  1.906252954  1.697092417  1.805815671  3.070084967  2.079309291  2.065501658  1.938043176
 [99]  1.190482099  1.095217145  3.495261507  1.785855790  1.576695253  1.685418507  2.949687803
[106]  1.958912127  1.945104495  1.817646012  1.070084935  0.974819981  3.596715295  1.887309579
[113]  1.678149042  1.786872296  3.051141592  2.060365916  2.046558284  1.919099801  1.171538724
[120]  1.076273770  3.514474264  1.805068548  1.595908011  1.704631265  2.968900561  1.978124885
[127]  1.964317253  1.836858770  1.089297693  0.994032739  3.699886880  1.990481164  1.781320627
[134]  1.890043881  3.154313177  2.163537501  2.149729868  2.022271386  1.274710309  1.179445355
[141]  3.453044373  1.743638657  1.534478120  1.643201374  2.907470670  1.916694994  1.902887361
[148]  1.775428879  1.027867802  0.932602848  3.500094790  1.790689074  1.581528537  1.690251791
[155]  2.954521087  1.963745411  1.949937778  1.822479296  1.074918219  0.979653265  1.662700613
[162] -0.046705103 -0.255865640 -0.147142386  1.117126910  0.126351234  0.112543601 -0.014914881
[169] -0.762475958 -0.857740912  1.618039603 -0.091366113 -0.300526650 -0.191803396  1.072465900
[176]  0.081690224  0.067882591 -0.059575891 -0.807136968 -0.902401922  1.497642439 -0.211763277
[183] -0.420923814 -0.312200560  0.952068736 -0.038706940 -0.052514572 -0.179973055 -0.927534132
[190] -1.022799086  1.599096228 -0.110309488 -0.319470025 -0.210746771  1.053522525  0.062746849
[197]  0.048939216 -0.078519266 -0.826080343 -0.921345297  1.516855197 -0.192550519 -0.401711056
[204] -0.292987802  0.971281494 -0.019494182 -0.033301815 -0.160760297 -0.908321374 -1.003586328
[211]  1.702267813 -0.007137903 -0.216298440 -0.107575186  1.156694110  0.165918434  0.152110801
[218]  0.024652319 -0.722908758 -0.818173712  1.455425306 -0.253980410 -0.463140947 -0.354417693
[225]  0.909851603 -0.080924073 -0.094731706 -0.222190188 -0.969751265 -1.065016219  1.502475723
[232] -0.206929993 -0.416090530 -0.307367276  0.956902020 -0.033873656 -0.047681289 -0.175139771
[239] -0.922700848 -1.017965802  1.218795228 -0.490610488 -0.699771026 -0.591047772  0.673221524
[246] -0.317554151 -0.331361784 -0.458820267 -1.206381343 -1.301646297  1.174134218 -0.535271498
[253] -0.744432036 -0.635708782  0.628560514 -0.362215161 -0.376022794 -0.503481277 -1.251042354
[260] -1.346307307  1.053737054 -0.655668662 -0.864829199 -0.756105945  0.508163351 -0.482612325
[267] -0.496419958 -0.623878440 -1.371439517 -1.466704471  1.155190843 -0.554214873 -0.763375411
[274] -0.654652156  0.609617140 -0.381158536 -0.394966169 -0.522424651 -1.269985728 -1.365250682
[281]  1.072949812 -0.636455904 -0.845616442 -0.736893188  0.527376109 -0.463399567 -0.477207200
[288] -0.604665682 -1.352226759 -1.447491713  1.258362428 -0.451043288 -0.660203826 -0.551480572
[295]  0.712788724 -0.277986951 -0.291794584 -0.419253067 -1.166814143 -1.262079097  1.011519921
[302] -0.697885795 -0.907046333 -0.798323079  0.465946217 -0.524829458 -0.538637091 -0.666095574
[309] -1.413656651 -1.508921604  1.058570338 -0.650835378 -0.859995916 -0.751272662  0.512996634
[316] -0.477779041 -0.491586674 -0.619045157 -1.366606234 -1.461871187

hist(predict(cox_fullmodel1))

Now we can compare the models to see which model is best, in this case it is fullmodel1.

anova(cox_simplemodel1,cox_fullmodel1) 

Analysis of Deviance Table
 Cox model: response is  Surv(NumDaysAlive, CensorReproduction)
 Model 1: ~Habitat * Treatment
 Model 2: ~Habitat * Treatment + (1 | Site) + (1 | Block)
   loglik  Chisq Df P(>|Chi|)    
1 -642.84                        
2 -623.87 37.938  2  5.78e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#See example: https://www.rdocumentation.org/packages/coxme/versions/2.2-16/topics/coxme
AIC(cox_simplemodel1,cox_fullmodel1) #But comparing the AIC scores is easiest. Keep the lower AIC score because that is considered the better model. Here it is the fullmodel1.

####Full Model 2 - coxme Now, let’s make a model with no interaction term and then we’ll compare that to the fullmodel1 which was deemed the best above.

cox_fullmodel2<-coxme(Surv(NumDaysAlive,CensorReproduction)~
                        Habitat+
                        Treatment+
                        (1|Site)+
                        (1|Block),
                      data=mydata) 
summary(cox_fullmodel2)

Cox mixed-effects model fit by maximum likelihood
  Data: mydata
  events, n = 160, 320
  Iterations= 8 53 
                    NULL Integrated    Fitted
Log-likelihood -676.4353  -624.8036 -607.6067

                   Chisq    df p    AIC   BIC
Integrated loglik 103.26  4.00 0  95.26 82.96
 Penalized loglik 137.66 12.28 0 113.10 75.34

Model:  Surv(NumDaysAlive, CensorReproduction) ~ Habitat + Treatment +      (1 | Site) + (1 | Block) 
Fixed coefficients
                                coef exp(coef)  se(coef)     z       p
HabitatVegetated         -0.02704818 0.9733143 0.1668362 -0.16 8.7e-01
TreatmentBiomass Removal  2.18264981 8.8697783 0.3276242  6.66 2.7e-11

Random effects
 Group Variable  Std Dev    Variance  
 Site  Intercept 0.15562304 0.02421853
 Block Intercept 0.81380754 0.66228271

Let’s compare the the first two models to test for the significance of the term that is removed (using LR). Here we see that fullmodel2 has a lower AIC score so now it is the best model.

anova(cox_fullmodel1,cox_fullmodel2) #Not significant

Analysis of Deviance Table
 Cox model: response is  Surv(NumDaysAlive, CensorReproduction)
 Model 1: ~Habitat * Treatment + (1 | Site) + (1 | Block)
 Model 2: ~Habitat + Treatment + (1 | Site) + (1 | Block)
   loglik  Chisq Df P(>|Chi|)
1 -623.87                    
2 -624.80 1.8627  1    0.1723

AIC(cox_fullmodel1,cox_fullmodel2) #The AIC scores are within 2 points of each other, so we can keep the simpler model, which is fullmodel2.

####Full Model 3 - coxme So, now let’s add in population nested under site as a random effect

cox_fullmodel3<-coxme(Surv(NumDaysAlive,CensorReproduction)~
                        Habitat+
                        Treatment+
                        (1|Site/Population)+
                        (1|Block),
                      data=mydata) 
summary(cox_fullmodel3)

Cox mixed-effects model fit by maximum likelihood
  Data: mydata
  events, n = 160, 320
  Iterations= 19 119 
                    NULL Integrated    Fitted
Log-likelihood -676.4353  -623.9201 -601.9647

                   Chisq    df p    AIC   BIC
Integrated loglik 105.03  5.00 0  95.03 79.65
 Penalized loglik 148.94 15.67 0 117.59 69.39

Model:  Surv(NumDaysAlive, CensorReproduction) ~ Habitat + Treatment +      (1 | Site/Population) + (1 | Block) 
Fixed coefficients
                                coef exp(coef)  se(coef)     z       p
HabitatVegetated         -0.01966282 0.9805292 0.2132455 -0.09 9.3e-01
TreatmentBiomass Removal  2.17096449 8.7667354 0.3310172  6.56 5.4e-11

Random effects
 Group           Variable    Std Dev     Variance   
 Site/Population (Intercept) 0.267339158 0.071470225
 Site            (Intercept) 0.019863810 0.000394571
 Block           Intercept   0.854306936 0.729840341

Now we can compare the models to see which model is best and we see that fullmodel3 has a lower AIC score.

anova(cox_fullmodel2,cox_fullmodel3) #Significant

Analysis of Deviance Table
 Cox model: response is  Surv(NumDaysAlive, CensorReproduction)
 Model 1: ~Habitat + Treatment + (1 | Site) + (1 | Block)
 Model 2: ~Habitat + Treatment + (1 | Site/Population) + (1 | Block)
   loglik  Chisq Df P(>|Chi|)
1 -624.80                    
2 -623.92 1.7671  1    0.1837

AIC(cox_fullmodel2,cox_fullmodel3) #Looks like fullmodel3 is the better model because of the lower AIC score

####Best Model The effect of source habitat was also not significant for the proportion of survival to bud (Z = -0.09, P = 0.93, Figure 4B). Again, the treatment effect was significant but not the interaction between treatment and source habitat (Z = 6.56, P < 0.0001).

summary(cox_fullmodel3)

Cox mixed-effects model fit by maximum likelihood
  Data: mydata
  events, n = 160, 320
  Iterations= 19 119 
                    NULL Integrated    Fitted
Log-likelihood -676.4353  -623.9201 -601.9647

                   Chisq    df p    AIC   BIC
Integrated loglik 105.03  5.00 0  95.03 79.65
 Penalized loglik 148.94 15.67 0 117.59 69.39

Model:  Surv(NumDaysAlive, CensorReproduction) ~ Habitat + Treatment +      (1 | Site/Population) + (1 | Block) 
Fixed coefficients
                                coef exp(coef)  se(coef)     z       p
HabitatVegetated         -0.01966282 0.9805292 0.2132455 -0.09 9.3e-01
TreatmentBiomass Removal  2.17096449 8.7667354 0.3310172  6.56 5.4e-11

Random effects
 Group           Variable    Std Dev     Variance   
 Site/Population (Intercept) 0.267339158 0.071470225
 Site            (Intercept) 0.019863810 0.000394571
 Block           Intercept   0.854306936 0.729840341

#####Risk Assessment These values are found in the model summary, but if you want to pull them out, here is how you interpret them. 1 = no effect, <1 = decreased risk of death, >1 = increased risk of death.

exp(coef(cox_fullmodel3)) #This should be interpreted that Biomass Removal is almost 9% more likely to survive to reproduction compared to Control.

        HabitatVegetated TreatmentBiomass Removal 
               0.9805292                8.7667354 

For the next manuscript: Looks like roadside and vegetated populations are the same, so let’s combine them together in a model (aka, removing the Habitat term)

cox_fullmodel4<-coxme(Surv(NumDaysAlive,CensorReproduction)~
                        Treatment+
                        (1|Site/Population)+
                        (1|Block),
                      data=mydata) 
summary(cox_fullmodel4)

Cox mixed-effects model fit by maximum likelihood
  Data: mydata
  events, n = 160, 320
  Iterations= 19 119 
                    NULL Integrated    Fitted
Log-likelihood -676.4353  -623.9245 -601.9447

                   Chisq    df p    AIC   BIC
Integrated loglik 105.02  4.00 0  97.02 84.72
 Penalized loglik 148.98 15.09 0 118.79 72.37

Model:  Surv(NumDaysAlive, CensorReproduction) ~ Treatment + (1 | Site/Population) +      (1 | Block) 
Fixed coefficients
                             coef exp(coef)  se(coef)    z       p
TreatmentBiomass Removal 2.172197  8.777544 0.3308417 6.57 5.2e-11

Random effects
 Group           Variable    Std Dev      Variance    
 Site/Population (Intercept) 0.2678093376 0.0717218413
 Site            (Intercept) 0.0198638048 0.0003945707
 Block           Intercept   0.8550652106 0.7311365144

summary(cox_fullmodel3)

Cox mixed-effects model fit by maximum likelihood
  Data: mydata
  events, n = 160, 320
  Iterations= 19 119 
                    NULL Integrated    Fitted
Log-likelihood -676.4353  -623.9201 -601.9647

                   Chisq    df p    AIC   BIC
Integrated loglik 105.03  5.00 0  95.03 79.65
 Penalized loglik 148.94 15.67 0 117.59 69.39

Model:  Surv(NumDaysAlive, CensorReproduction) ~ Habitat + Treatment +      (1 | Site/Population) + (1 | Block) 
Fixed coefficients
                                coef exp(coef)  se(coef)     z       p
HabitatVegetated         -0.01966282 0.9805292 0.2132455 -0.09 9.3e-01
TreatmentBiomass Removal  2.17096449 8.7667354 0.3310172  6.56 5.4e-11

Random effects
 Group           Variable    Std Dev     Variance   
 Site/Population (Intercept) 0.267339158 0.071470225
 Site            (Intercept) 0.019863810 0.000394571
 Block           Intercept   0.854306936 0.729840341

anova(cox_fullmodel3,cox_fullmodel4) #Significant

Analysis of Deviance Table
 Cox model: response is  Surv(NumDaysAlive, CensorReproduction)
 Model 1: ~Habitat + Treatment + (1 | Site/Population) + (1 | Block)
 Model 2: ~Treatment + (1 | Site/Population) + (1 | Block)
   loglik  Chisq Df P(>|Chi|)
1 -623.92                    
2 -623.92 0.0088  1    0.9253

AIC(cox_fullmodel3,cox_fullmodel4) #Using the full dataset (not loaded here) it looks like fullmodel4 is the better model because the AIC score is within 2 points of each other, therefore the models are assessed the same and you should take the simpler model. But this is for another manuscript :-)

#Reproductive Biomass This code uses Biomass data with Habitat (roadside and vegetated) and Treatment in a lmer model. ANOVA and Tukey are used on the successful Model 3 with the creation of a box plot as a finished product.

Note: 10 blocks, 5 treatments, 16 populations (CHE-O), 8 sites (random: population pairs; CHE), 2 habitat (fixed effect: roadside and vegetated; R or O)

Also Note: Biomass data is a measure of aboveground biomass of all individuals harvested from the site. In some cases this was before they budded, but we harvested the wilted plant in case that information was needed in the future. Here we only want to look at biomass (well, log(Biomass)) for reproductive individuals.

Here we need to subset the data to only look at Biomass when CensorReproduction = 1, so that all Biomass data is for reproductive individuals only.

reproduction_mydata<-subset(mydata,CensorReproduction%in%c('1')) #Here I am only looking at the Biomass data where the CensorReproduction = 1
repro_mydata<-na.omit(reproduction_mydata)

##Histograms ###Original data The original data is skewed and fails the Shapiro-Wilk normality test, so we should consider a log transformation.

#All data
hist(repro_mydata$Biomass,col='steelblue',
     main='Original') #Original data is skewed, let's test for normality and consider log transforming the data

shapiro.test(repro_mydata$Biomass) #fails the Shapiro-Wilk normality test

    Shapiro-Wilk normality test

data:  repro_mydata$Biomass
W = 0.7266, p-value = 9.233e-16

#Biomass Removal data
biomass_hist<-repro_mydata %>%
  select(Biomass,Treatment) %>% 
  filter(Treatment =="Biomass Removal")
hist(biomass_hist$Biomass,breaks=20,col='steelblue',main='Original Biomass Removal')

###Log transform data (https://www.statology.org/transform-data-in-r/) Log transformed data has a better distribution than the original data (although it still fails the Shapiro-Wilk normality test) so we will use the log transformed data with our models.

log_reproduction_mydata<-log10(repro_mydata$Biomass)
hist(log_reproduction_mydata,col='steelblue',
     main='Log Transformed') #Log transformed data, this looks better than the original distribution

shapiro.test(log_reproduction_mydata) #but it fails this test

    Shapiro-Wilk normality test

data:  log_reproduction_mydata
W = 0.94952, p-value = 1.959e-05

qqnorm(log_reproduction_mydata) #qqplot
qqline(log_reproduction_mydata) #add the line... kinda wobbly around the ends

###Square root transform data Nope, not square root

sqrt_biomass_mydata<-sqrt(repro_mydata$Biomass)
hist(sqrt_biomass_mydata,
     col='steelblue',main='Square Root Transformed') #nope, this does not look better

shapiro.test(sqrt_biomass_mydata)

    Shapiro-Wilk normality test

data:  sqrt_biomass_mydata
W = 0.93708, p-value = 1.987e-06

qqnorm(sqrt_biomass_mydata) #qqplot
qqline(sqrt_biomass_mydata) #add the line 

###Poisson distribution Poisson is actually not a good fit because the growth rate is not an integer (has decimals), which poisson cannot deal with. Maybe I’ll try gamma?

library(MASS)
MASS::fitdistr((repro_mydata$Biomass*1000),
               "Poisson")

     lambda   
  6420.076433 
 (   6.394701)

qqPlot(repro_mydata$Biomass,
       distribution="pois",
       lambda=1)

[1]  57 108

###Gamma distribution Maybe?

gamma<-fitdistr(repro_mydata$Biomass,
                "gamma")
qqp(repro_mydata$Biomass,
    "gamma",
    shape=gamma$estimate[[1]],
    rate=gamma$estimate[[2]])

[1]  57 108

##Models Julia recommends the log and gamma… then use DHARMa to test the dispersion of each models

###Full Model 1 - lmer: Modeling Habitat and Treatment with Site and Block as random effects Site as a random effect does not explain any of the variance in the model, therefore let’s try Site as a fixed effect to demonstrate that it doesn’t add to the model.

fullmodel1<-lmer(log(Biomass)~
                   Habitat*Treatment+
                   #(1|Site)+ we dropped site because it was a singular fit because it is described by habitat
                   (1|Block),
                 data=repro_mydata)
isSingular(fullmodel1,tol=1e-4) #false without site

[1] FALSE

summary(fullmodel1) #but this runs fine

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log(Biomass) ~ Habitat * Treatment + (1 | Block)
   Data: repro_mydata

REML criterion at convergence: 414.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2800 -0.5732  0.0033  0.6569  2.4381 

Random effects:
 Groups   Name        Variance Std.Dev.
 Block    (Intercept) 0.2717   0.5212  
 Residual             0.7152   0.8457  
Number of obs: 157, groups:  Block, 10

Fixed effects:
                                          Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                                -1.4978     0.2579  35.5278  -5.808 1.31e-06 ***
HabitatVegetated                            0.5960     0.2915 145.3124   2.045   0.0427 *  
TreatmentBiomass Removal                    3.1903     0.2271 145.9753  14.049  < 2e-16 ***
HabitatVegetated:TreatmentBiomass Removal  -0.6511     0.3284 144.7802  -1.982   0.0493 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) HbttVg TrtmBR
HabittVgttd -0.512              
TrtmntBmssR -0.675  0.578       
HbttVgt:TBR  0.455 -0.884 -0.671

anova(fullmodel1)

Type III Analysis of Variance Table with Satterthwaite's method
                   Sum Sq Mean Sq NumDF  DenDF  F value  Pr(>F)    
Habitat             1.913   1.913     1 145.42   2.6746 0.10412    
Treatment         206.045 206.045     1 147.25 288.0910 < 2e-16 ***
Habitat:Treatment   2.811   2.811     1 144.78   3.9301 0.04932 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

predict(fullmodel1)

          1           4           5           6           7           8           9          10 
 0.95427788  2.20663461  1.32798895  2.05847024  1.46081666  1.84178239  1.65887219  1.20629641 
         11          13          14          15          16          18          19          21 
 0.95427788  1.68228142  2.20663461  1.32798895  2.05847024  1.84178239  1.65887219  0.95427788 
         22          23          24          25          26          27          28          29 
 2.52743203  1.68228142  2.20663461  1.32798895  2.05847024  1.46081666  1.84178239  1.65887219 
         30          31          32          33          35          37          38          40 
 1.20629641  0.95427788  2.52743203  1.68228142  1.32798895  1.46081666  1.84178239  1.20629641 
         41          42          44          46          47          48          49          51 
 0.95427788  2.52743203  2.20663461  2.05847024  1.46081666  1.84178239  1.65887219  0.95427788 
         52          53          55          56          57          58          59          60 
 2.52743203  1.68228142  1.32798895  2.05847024  1.46081666  1.84178239  1.65887219  1.20629641 
         62          63          64          66          67          68          69          71 
 2.52743203  1.68228142  2.20663461  2.05847024  1.46081666  1.84178239  1.65887219  0.95427788 
         72          73          76          77          78          81          83          85 
 2.52743203  1.68228142  2.05847024  1.46081666  1.84178239  0.89922077  1.62722431  1.27293184 
         86          87          88          89          90          93          94          95 
 2.00341312  1.40575955  1.78672528  1.60381508  1.15123929  1.62722431  2.15157750  1.27293184 
         96          97          98         100         102         103         105         106 
 2.00341312  1.40575955  1.78672528  1.15123929  2.47237491  1.62722431  1.27293184  2.00341312 
        107         109         111         112         113         114         116         118 
 1.40575955  1.60381508  0.89922077  2.47237491  1.62722431  2.15157750  2.00341312  1.78672528 
        119         120         122         123         124         127         128         129 
 1.60381508  1.15123929  2.47237491  1.62722431  2.15157750  1.40575955  1.78672528  1.60381508 
        130         131         133         134         135         136         137         138 
 1.15123929  0.89922077  1.62722431  2.15157750  1.27293184  2.00341312  1.40575955  1.78672528 
        139         140         141         142         143         144         147         148 
 1.60381508  1.15123929  0.89922077  2.47237491  1.62722431  2.15157750  1.40575955  1.78672528 
        149         151         153         154         155         156         157         158 
 1.60381508  0.89922077  1.62722431  2.15157750  1.27293184  2.00341312  1.40575955  1.78672528 
        159         160         161         163         164         170         172         174 
 1.60381508  1.15123929 -2.23601152 -1.50800798 -0.98365479 -1.98399300 -0.66285738 -0.98365479 
        184         185         189         195         196         201         210         211 
-0.98365479 -1.86230045 -1.53141721 -1.86230045 -1.13181917 -2.23601152 -1.98399300 -2.23601152 
        215         224         226         228         236         254         259         260 
-1.86230045 -0.98365479 -1.13181917 -1.34850702 -1.13181917 -0.38761839 -0.93538081 -1.38795660 
        270         272         280         281         284         285         286         287 
-1.38795660 -0.06682098 -1.38795660 -1.63997513 -0.38761839 -1.26626406 -0.53578277 -1.13343635 
        295         296         300         303         310 
-1.26626406 -0.53578277 -1.38795660 -0.91197159 -1.38795660 

hist(predict(fullmodel1,type="response"))

qqnorm(resid(fullmodel1)) #qqplot
qqline(resid(fullmodel1)) #add the line

testDispersion(fullmodel1) #red line should be in the middle of the distribution

    DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated

data:  simulationOutput
dispersion = 0.98467, p-value = 0.968
alternative hypothesis: two.sided

myDHARMagraph1<-simulateResiduals(fullmodel1) #making a graph using DHARMa package, also testing for heteroscedasticity, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html#heteroscedasticity
plot(myDHARMagraph1) #plotting graph. At this point, you don't want any text or lines to be red.

###Full Model 2 - glmm

fullmodel2<-glmmTMB(Biomass~
                             Habitat*Treatment+
                             (1|Site)+
                             (1|Block),
                           family=Gamma(link="log"),
                           data=repro_mydata) 
summary(fullmodel2)

 Family: Gamma  ( log )
Formula:          Biomass ~ Habitat * Treatment + (1 | Site) + (1 | Block)
Data: repro_mydata

     AIC      BIC   logLik deviance df.resid 
   738.8    760.2   -362.4    724.8      150 

Random effects:

Conditional model:
 Groups Name        Variance  Std.Dev. 
 Site   (Intercept) 1.505e-09 3.879e-05
 Block  (Intercept) 1.794e-01 4.236e-01
Number of obs: 157, groups:  Site, 8; Block, 10

Dispersion estimate for Gamma family (sigma^2): 0.583 

Conditional model:
                                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)                                -1.1677     0.2258  -5.171 2.32e-07 ***
HabitatVegetated                            0.6390     0.2636   2.424  0.01536 *  
TreatmentBiomass Removal                    3.2082     0.2128  15.074  < 2e-16 ***
HabitatVegetated:TreatmentBiomass Removal  -0.7766     0.2975  -2.610  0.00904 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

qqnorm(resid(fullmodel2)) #qqplot
qqline(resid(fullmodel2)) #add the line

testDispersion(fullmodel2) #red line should be in the middle of the distribution

    DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated

data:  simulationOutput
dispersion = 1.2635, p-value = 0.48
alternative hypothesis: two.sided

myDHARMagraph2<-simulateResiduals(fullmodel2) #making a graph using DHARMa package, also testing for heteroscedasticity, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html#heteroscedasticity
plot(myDHARMagraph2) #plotting graph. Looks good, but check the outliers to make sure they are real.

###Best Model

fullmodel1<-lmer(log(Biomass)~
                   Habitat*Treatment+
                   #(1|Site)+ we dropped site because it was a singular fit because it is described by habitat
                   (1|Block),
                 data=repro_mydata)
fullmodel1

Linear mixed model fit by REML ['lmerModLmerTest']
Formula: log(Biomass) ~ Habitat * Treatment + (1 | Block)
   Data: repro_mydata
REML criterion at convergence: 414.1645
Random effects:
 Groups   Name        Std.Dev.
 Block    (Intercept) 0.5212  
 Residual             0.8457  
Number of obs: 157, groups:  Block, 10
Fixed Effects:
                              (Intercept)                           HabitatVegetated  
                                  -1.4978                                     0.5960  
                 TreatmentBiomass Removal  HabitatVegetated:TreatmentBiomass Removal  
                                   3.1903                                    -0.6511  

#ggplot - interaction plots Treatment (Biomass Removal and Grassland) on the x-axis, using Roadside and Vegetated data

##Growth Rate x Treatment: Growth

pd<-position_dodge(0)

growth_mydata2<-growth_mydata%>% 
  replace_na(list(Growth.Rate=0))%>%
  filter(CensorReproduction==1)%>%
  #filter(Treatment=="Biomass Removal"|Treatment=="Grassland")%>%
  group_by(Treatment,Habitat)%>%
  summarise(Mean=mean(Growth.Rate),
            SD=sd(Growth.Rate),
            N=length(Growth.Rate))%>%
  mutate(SE=SD/sqrt(N))

`summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.

field_int_rate_grow_gg<-ggplot(growth_mydata2,
                        aes(x=factor(Treatment,levels=c("Biomass Removal","Grassland")),
                            y=Mean,
                            shape=Habitat,
                            group=Habitat,
                            color=Habitat))+
  coord_cartesian(ylim = c(0,0.5))+
  xlab("")+
  theme_bw(base_size=16)+
  theme(panel.border=element_blank(),
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        axis.line=element_line(colour="black"),
        axis.title.y=ggtext::element_markdown(),
        legend.position = c(0.8, 0.8))+
 # ggtitle("Growth of Dittrichia graveolens that budded")+
  ylab("*Dittrichia graveolens* \n (Change in Leaf Length/Days)")+
  geom_line(size=0.8,
            position=pd)+
  geom_point(size=4,
             position=pd)+
  scale_shape_manual(values=c(16,2))+
  scale_color_manual(values=c("gray85","forestgreen"))+
  geom_errorbar(width=0.15,
                position=pd,
                aes(ymin=(Mean-SE),
                    ymax=(Mean+SE)))+
  labs(y="Change in Leaf Length (mm/day)",
       color="Source Habitat",
       shape="Source Habitat")
field_int_rate_grow_gg

#ggsave(plot=field_int_rate_grow_gg,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/field_int_rate_grow_gg.png",width=25,height=15,units="cm",dpi=800)

##Survival to bud X Treatment

pd<-position_dodge(0)

surv_mydata1<-mydata%>%
  #filter(Treatment=="Biomass Removal"|Treatment=="Grassland")%>%
  group_by(Treatment,Habitat)%>%
  summarise(Mean=mean(PropBudSite),
            SD=sd(PropBudSite),
            N=length(PropBudSite))%>%
  mutate(SE=SD/sqrt(N))

`summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.

field_int_surv_all_gg<-ggplot(data=surv_mydata1,
                        aes(x=factor(Treatment,levels=c("Biomass Removal","Grassland")),
                            y=Mean,
                            shape=Habitat,
                            group=Habitat,
                            color=Habitat))+
  ylab("Proportion Survival to Bud")+
  coord_cartesian(ylim=c(0,1))+
  xlab("")+
  theme_classic(base_size=16)+
  theme(panel.border=element_blank(),
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        axis.line=element_line(colour="black"),
        legend.position = c(0.8, 0.8))+
  scale_y_continuous(name="Proportion Survival to Bud",
                     limits=c(0,1),
                     breaks=c(0.0,0.2,0.4,0.6,0.8,1.0))+
  geom_line(size=0.8,
            position=pd)+
  geom_point(size=4,
             position=pd)+
  scale_shape_manual(values=c(16,2))+
  scale_color_manual(values=c("gray85","forestgreen"))+
  geom_errorbar(width=0.15,
                position=pd,
                aes(ymin=(Mean-SE),
                    ymax=(Mean+SE)))+
  labs(y="Proportion Survival to Bud",
       fill="Source Habitat",
       color="Source Habitat",
       shape="Source Habitat")
field_int_surv_all_gg

#ggsave(plot=field_int_surv_all_gg,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/field_int_surv_all_gg.png",width=25,height=15,units="cm",dpi=800)

##Biomass x Treatment: Growth

pd<-position_dodge(0)

growth_mydata1<-growth_mydata%>% 
  replace_na(list(Biomass=0))%>%
  filter(CensorReproduction==1)%>%
#  filter(Treatment=="Biomass Removal"|Treatment=="Grassland")%>%
  group_by(Treatment,Habitat)%>%
  summarise(Mean=mean(Biomass),
            SD=sd(Biomass),
            N=length(Biomass))%>%
  mutate(SE=SD/sqrt(N))

`summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.

field_int_bio_grow_gg<-ggplot(growth_mydata1,
                        aes(x=factor(Treatment,levels=c("Biomass Removal","Grassland")),
                            y=Mean,
                            shape=Habitat,
                            group=Habitat,
                            color=Habitat))+
  ylab("Biomass (g)")+
  coord_cartesian(ylim = c(0,10))+
  xlab("")+
  #ggtitle("Biomass of Dittrichia graveolens that budded")+
  theme_bw(base_size=16)+
  theme(panel.border=element_blank(),
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        axis.line=element_line(colour="black"),
        axis.title.y=ggtext::element_markdown(),
        legend.position = c(0.8, 0.8))+
  geom_line(size=0.8,
            position=pd)+
  geom_point(size=4,
             position=pd)+
  scale_shape_manual(values=c(16,2))+
  scale_color_manual(values=c("gray85","forestgreen"))+
  geom_errorbar(width=0.15,
                position=pd,
                aes(ymin=(Mean-SE),
                    ymax=(Mean+SE)))+
  labs(y="*D. graveolens* Biomass (g)",
       color="Source Habitat",
       shape="Source Habitat")
field_int_bio_grow_gg

#ggsave(plot=field_int_bio_grow_gg,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/field_int_bio_grow_gg.png",width=25,height=15,units="cm",dpi=800)

##Figure 4 panel (3 plots)

library(patchwork)

multi_panel_gg<-(field_int_rate_grow_gg+theme(legend.position = c(0.8,0.9)))/
  (field_int_surv_all_gg+theme(legend.position="none"))/
  (field_int_bio_grow_gg+theme(legend.position="none"))+
  plot_annotation(tag_levels=c("A"),
                  tag_suffix=")  ")&
  theme(text=element_text(size=8))

multi_panel_gg

ggsave(plot=multi_panel_gg,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/multi_panel_gg.png",width=8.5,height=17,units="cm",dpi=800)

#Survival to Bud ##Models ###Best Model

exp(coef(cox_fullmodel4))

                         HabitatVegetated                  TreatmentBiomass Removal HabitatVegetated:TreatmentBiomass Removal 
                                0.9850102                                 0.1762327                                 0.9785911 

#ggplot - survival curves ##Figure 5

###autoplot

#Now let's try making survival curves to plot, however, this is not the right model for this figure
cox_model_graph3<-survfit(Surv(NumDaysAlive,CensorReproduction)~
                            Habitat+
                            Treatment,
                          data=mydata)

surv_habitat_plot1<-autoplot(cox_model_graph3)+
  labs(x="Survival (Days)",
       y="Probability")+
  theme(plot.title=element_text(hjust=0.5),
        axis.title.x=element_text(face="bold",
                                  color="Black",
                                  size=12),
        axis.title.y=element_text(face="bold",
                                  color="Black",
                                  size=12))
surv_habitat_plot1
ggsave(plot=surv_habitat_plot1,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/surv_habitat_plot1.png",width=16,height=8,units="cm",dpi=800)

cox_model_graph4<-survfit(Surv(NumDaysAlive,CensorReproduction)~
                            Treatment+
                            Habitat,
                          data=mydata)

surv_habitat_plot2<-autoplot(cox_model_graph4)+
  labs(x="Survival (Days)",
       y="Probability")+
  theme(plot.title=element_text(hjust=0.5),
        axis.title.x=element_text(face="bold",
                                  color="Black",
                                  size=12),
        axis.title.y=element_text(face="bold",
                                  color="Black",
                                  size=12))
surv_habitat_plot2
ggsave(plot=surv_habitat_plot2,file="/Users/Miranda/Documents/Education/UC Santa Cruz/Dittrichia/Dittrichia_Analysis/figures/surv_habitat_plot2.png",width=16,height=8,units="cm",dpi=800)

New plot idea http://sthda.com/english/wiki/survminer-r-package-survival-data-analysis-and-visualization#change-legend-title-labels-and-position

https://www.rdocumentation.org/packages/survminer/versions/0.4.9#computing-and-passin-p-values

###ggsurvplot